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:3 ABSTRACT 

^^ We carried out extensive numerical orbit integrations to probe the long-term chaotic 

,-C| dynamics of the two strongest mean motion resonances of Neptune in the Kuiper belt, 

the 3:2 (Plutinos) and 2:1 (Twotinos). Our primary results include a computation of the 
relative volumes of phase space characterized by large- and small-resonance libration 



Ph 



o 



t/j amplitudes, and maps of resonance stability measured by mean chaotic diffusion rate. 

I ^1 We find that Neptune's 2:1 resonance has weaker overall long-term stability than the 

3:2 — only ~ 15% of Twotinos are projected to survive for 4 Gyr, compared to ^ 27% 

j>. of Plutinos, based on an extrapolation from our 1-Gyr integrations. We find that 

^*^ Pluto has only a modest effect, causing a ~ 4% decrease in the Plutino population 

QQ that survives to 4 Gyr. Given current observational estimates, and assuming an initial 

^N distribution of particles proportional to the local phase space volume in the resonance, 

^^ we conclude that the primordial populations of Plutinos and Twotinos formerly made 

QQ up more than half the population of the classical and resonant Kuiper Belt. We also 

^^ conclude that Twotinos were originally nearly as numerous as Plutinos; this is consistent 

J>" with predictions from early models of smooth giant planet migration and resonance 

S^ sweeping of the Kuiper Belt, and provides a useful constraint for more detailed models. 

Subject headings: Celestial Mechanics, Kuiper Belt 



1. Introduction 

Over the past 15 years of observations, it has become clear that the trans-Neptune population 
of Kuiper Belt objects (KBOs) consists of several distinct dynamical classes ( Jewitt and Luu||2000 



Malhotraet al.|2000[ Gladman et al.|2008 1 . The "classical" KBOs, found primarily in the semimajor 
axis range of 40 AU to 48 AU, consist of two sub-classes: a dynamically cold population having 
relatively circular and low-inclination orbits, and a dynamically hot population with higher orbital 



eccentricities and inclinations ( Brown|2001 Levison and Stern|2001 Trujillo and Brown|2002 1. The 
scattered disk objects (SDOs) have highly eccentric and inclined orbits, with perihelia a few AU 
beyond Neptune's orbit ( [Duncan and Levison 1997); a small number of the SDOs with perihelia 
> 40 AU are thought to be a distinct population, the "extended scattered disk" (also known as 



the 'detached objects') (Gladman et al. 2002 1. Finally, there are the resonant KBOs, which have 
orbits in mean motion resonance with Neptune; approximately 23% of all known KBOs are in 



this class (Chiang et al. 2007). Most of the known resonant KBOs are in the 3:2 resonance at 



semimajor axis a = 39.4 AU, and have been dubbed 'Plutinos' in recognition of the largest and 
longest-observed member of this population; the 2:1 resonance at a = 47.7 AU has the next-highest 
observed resonant population, and have been dubbed 'Twotinos'; several other resonances are also 



populated (Chiang et al. 2003 2007). 



The particular interest of this paper is the dynamics of bodies within the 3:2 and 2:1 mean- 
motion resonances with Neptune. Because the structure of these resonances provides a readily 
understandable mechanism for preserving a population for 4 Gyr or more, as well as dynamical 
pathways to less stable orbits where they may encounter Neptune and be transported elsewhere in 
the solar system, resonant KBOs have been proposed as the source populations for the Centaurs and 
Jupiter- family comets ( Duncan et al.|[l995 Duncan and Levison||1997 Morbidelli 1997). Resonant 
KBOs are also of great interest because their relative populations hold clues to the orbital migration 



history of the giant planets and the dynamical history of the outer solar system as a whole (Malhotra 



1995 MorbideUi and Brown 2004; Morbidelh et al. 2008). 



In the present work, we use numerical analysis to obtain a detailed map of stability as it varies 
with location in the phase space of the 3:2 and 2:1 mean- motion resonances with Neptune. One 
major goal of this study is to estimate (by extrapolating the rates seen in our 1-Gyr numerical 
integrations) the Plutino and Twotino populations ~ 4 Gyr ago when the dynamical structure of 
the Kuiper belt was presumably established. We also investigate the effects of Pluto on the Plutino 
population, and the behavior of particles after they have escaped from resonance. 

Section [2] of this paper discusses the properties of resonant particles and our treatment of 
chaotic diffusion; Section [3] describes the numerical experiments that we carried out to study the 
long-term stability of particles in the 3:2 and 2:1 resonances; Section [4] describes the results of our 
numerical models; Section [5] provides discussion, analysis, and comparison to previous work; and 
Section [6] gives a summary and conclusion. 



2. Resonances and Chaotic Diffusion 



The dynamics of resonant KBOs are characterized by the libration amplitude of their resonant 
argument, (j). For the 3:2 and the 2:1 exterior mean motion resonances of Neptune, we will be 



concerned with the particular resonant arguments, 



^3;2 = 3A — 2AAr — w, 
^2:1 = 2A — Xn — W, 



:i) 



where Atv is the mean longitude of Neptune, and A and zu are the mean longitude and longitude of 
perihelion, respectively, of a KBO. 

There is one important difference between the 3:2 and the 2:1 resonance librations. The reso- 
nant arguments of particles in most exterior resonances librate about 180°; however, the libration 
center of a l:n resonant argument is eccentricity-dependent and bifurcated into multiple asymmet- 



ric libration centers (Beauge 1994 Malhotra 1996 Winter and Murray 1997; Pan and Sari 2004 1 . 



For the 2:1 resonance in particular, the asymmetric libration centers appear as a bifurcation of the 
exact resonant circular periodic orbit at (/> = 180°; for e ~ 0.1, the libration centers are near ±90°, 
and asymptotically approach values near (f> = ±60° for increasing values of eccentricity (see figure 
in Malhotra 19991. A particle whose resonant argument (j) librates about the positive value (0c) 
reaches its perihelion near a longitude that is 0c ahead of Neptune's position at that time; thus such 
particles are referred to as "leading" librators. Conversely, a particle whose (p librates about the 
corresponding negative value (—0c) reaches its perihelion at a longitude that is 0c behind Neptune's 
position, and is thus referred to as a "trailing" librator. "Symmetric" librators, whose librates 
about 180°, can continue to exist even when the asymmetric libration centers are present, but they 
necessarily have large libration amplitudes because they must effectively transition from the leading 



lobe to the trailing lobe and back to the leading lobe during each libration cycle (Malhotra 19961; 
these are analogous to "horseshoe orbits" in the classical three-body problem. 

The amplitude of the resonance libration is an important indicator of the stability of a reso- 
nant orbit, because a close encounter with Neptune becomes more likely as the particle's resonant 
argument strays away from the center of libration. We measure this libration amplitude, A0, as 
half the difference between the minimum and maximum excursions of 0. In the simplest analytical 
models for resonance, the libration amplitude is a constant of the motion. However, in reality, A0 is 
not an exact constant because additional weak resonances and near-resonances, associated not only 
with Neptune but with the other planets (primarily Jupiter, Saturn, and Uranus), interact with the 
dominant resonant perturbation from Neptune and produce quasi-periodic perturbations as well as 



chaotic behavior (e.g. Lecar et al. 2001 1 . To assess the long term stability of resonant KBOs, we 
are particularly interested in the chaotic evolution that may lead to escape from resonance. 

There are two defining characteristics of a chaotic dynamical system. One of these is the 
local exponential divergence of initially nearby particle orbits, and is quantified by the Lyapunov 



exponent (e.g., see Nesvorny and Roig 2000); however, the relationship between the Lyapunov 
exponent and escape time is not a simple one. The second characteristic of dynamical chaos is 
the existence of broad frequency bands in the Fourier spectrum of the time series of a dynamical 
variable (in contrast with the line spectrum for quasi-periodic regular dynamics). This property 
has been used to detect the presence of chaos and to measure the sizes of chaotic zones in the 



dynamics of the major planets of the solar system (e.g. Laskar 1990 1 . However, this also does not 
provide a direct estimate of the escape times that are of interest here. For the purpose of the 
present paper, a more practical approach approximates the chaotic behavior as a classical diffusion 
process as follows. We can define the approximate constants of the motion - the proper elements 
- over a time interval. At, which is long compared to the quasiperiodic perturbations but small 
compared to the long time of interest. (In a low order analytical perturbation theory, these proper 
elements are exact constants of the motion, but their constancy is not guaranteed when higher order 
perturbations are included.) The changes of the proper elements over successive time intervals can 
be attributed to chaos when one measures an increase with time of the dispersion of the proper 
elements; this is referred to as "chaotic diffusion". This approach has been developed by ILaskar] 



(1994) and was adopted by Morbidelli (19971 to study the chaotic evolution of Plutinos. We adopt 
this approach in the present work, as it lends itself to a clear analysis and interpretation of the long 
term changes in the populations of resonant KBOs. 



Following Morbidelli (19971, we define proper elements that are derived directly from the oscu- 
lating orbital elements by means of a numerical procedure described in Section |3.2[ this procedure 
removes quasi-periodic variations so that the chaotic behavior can be isolated. As in [Morbidelir 



( 1997 1 's study of the Plutinos, we observe 4 classes of behavior among the particles in our numerical 
models, examples of which are shown in Fig. [T| A particle's proper semimajor axis may remain 
nearly stationary over the integration, either because the particle is on a stable trajectory (an "in- 
variant torus" in phase space) or because its chaotic evolution is below the numerical resolution of 
our model (Fig. flk); secondly, it may evolve over a significant distance in phase space but remain 
in resonance (Fig. IIV)); thirdly, it may evolve within the resonance far enough to reach the chaotic 
zone and thence escape to a non-resonant orbit (Fig. lip) ; or fourthly, it may be strongly chaotic 
from the start and escape the resonance almost immediately (Fig. ^). 

In order to quantify the chaotic behavior for initial conditions across the resonance phase space, 
we assume that the chaotic evolution of the proper elements can be approximated as a diffusive 
process. In a diffusive process, ensembles of particles that start out nearby in phase space drift 
apart, as their proper elements execute a random walk. The rate at which they disperse is described 
by a diffusion coefficient, which we define as D = ((Aa)^)/At, where Aa is the change in proper 
semimajor axis over a time interval At, and the mean is taken over the ensemble of particles located 
in a small volume of phase space. The diffusion coefficient may vary with location in phase space, 
in which case particles will tend to spend more time in regions with small diffusion coefficient, and 



less time in regions with large diffusion coefficient. As explained in Section 3.2 below, we choose a 
"running window" time interval At = 10 Myr, which is much longer than any resonance libration 
periods, but much shorter than the ~ 10^ Myr timescale of interest for the history of the resonant 
KBOs. 
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Fig. 1. — The evolution of four selected particles from our models illustrates the classes of behavior 
that we observe: a) stable resonant orbit with little measurable diffusion, b) stable resonant orbit 
which diffuses within resonance, c) initially resonant orbit which eventually diffuses out of resonance 
(after 809 Myr in this case), d) strongly chaotic trajectory which quickly leaves resonance (after 



37 Myr in this case). The minimum proper elements (see Section 3.2) are plotted here, so the 
right-hand side of each plot is the center of the 3:2 resonance. Arrows point to the initial values. 
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3. Numerical Model and Analysis 

3.1. Pre-runs and Initial Conditions 

To obtain a sample of initial conditions that comprehensively covers the resonance zone, a 
series of "pre-runs" was performed. The primary purpose of the pre-runs is to identify values of the 
initial conditions that result in resonant particles, and then to include only those particles in the 
full 1-Gyr integration. This step is necessary because the initial osculating orbital elements do not 
have a simple relationship with the proper elements that define the parameters of the resonance 
phase space; this complexity owes to the large-amplitude short-period variations caused by the 
giant planets (Section |3.2[ ). For the full 1-Gyr runs, we then chose a subset of initial conditions 
from the pre-runs to ensure adequate coverage of the entire resonance zone, including the deepest 
zone at very small libration amplitudes; the latter has usually been poorly covered in previous 
studies. 

The orbital intergrations were performed using the "Swift-Skeel" mixed-variable symplectic 



A^-body integrator (Duncan et al. 1998 Wisdom and Holman 1991); we used a step size of 0.5 



Earth years. All runs included the 4 giant planets, whose initial orbital elements (obtained from 
JPlFj) corresponded to a starting epoch of 1997 June 1. The mutual perturbations of the massive 
planets were fully accounted for. 

For the pre-runs, our set of test particles began with 20 eccentricity values, evenly spaced in 
the interval e = [.05, .5]; and 7 inclination values Jm evenly spaced in the interval i = [5°, 35°]. For 
Plutinos, an 8th inclination value of 17.5° was added to increase the resolution in Pluto's vicinity. 
Each {e, i} bin contained 250 particles, with initial values of the angular elements {to, 0,, and M) 
randomly distributed in the interval M = [0°,360°). The total number of particles was 35,000 
Twotinos, and 40,000 Plutinos. All particles began with semimajor axis a = 39.35 AU for Plutinos, 
a = 47.87 AU for Twotinos. The orbits of these particles were then integrated, along with the 
giant planets, for 150,000 yr; this is several times the typical resonance libration period, and long 
enough for stable libration of the resonant argument to manifest itself. At each 2,000-yr interval, 
we recorded for each planet and test particle the minimum and maximum values of a, e, and i over 
that interval, as well as instantaneous values of 0,, uj, and M. 

For each particle in the pre-run, we track the evolution of the resonant argument, 4>s:2 for 
Plutinos and (j)2:i for Twotinos (see Eq. Q; and we verify persistence of libration by visual exam- 
ination of the time series plots of these resonant arguments. For librating particles, we calculate 
the libration amplitude of the resonant argument, A(j) = {4>max — 4'min)/'^- We found that the 
75 instantaneous points over the 150-kyr interval were sufficient to characterize the libration. For 



* http://ssd.jpl. nasa.gov/horizons. html 

'In this paper, we measured inclinatioi 
the Laplace plane of the Kuiper Belt (i.e., the plane about which orbital planes actually precess) is 1.86° ([Brown and Pan] 



'In this paper, we measured inclinations with respect to the J2000 ecliptic plane. The angle between the ecliptic plane and 



20041, and the mean correction is < 1° (Trujillo and Brown|2002l, much smaller than the inclinations discussed in this paper. 



Twotinos, we additionally categorized each particle, by visual inspection of the time series of (j), as 
having stable libration in one of three modes: leading ((/> < 180° at all times), trailing (0 > 180° 
at all times), or symmetric (0 librates about 180°). Particles that did not exclusively exhibit one 
of these behaviors during the pre-run were discarded as they are of no interest for the long term 
history of the resonant KBOs. 

Initial particles for the full 1-Gyr integration were then chosen from among the particles that 
exhibited stable libration in the pre-run. One in 16 particles with A(f> > 30° was randomly selected. 
For particles with A(f) < 30°, the probability of selection was increased to one in 4, to enhance the 
resolution of our results for tightly-bound librators; these particles were weighted by 0.25 in the 
subsequent analysis. The result was a set of 1,331 Plutinos and 1,445 Twotinos in the full 1-Gyr 
integration. The selected Twotinos are broken down as 479 leading, 636 trailing, and 330 symmetric 
librators; the weighted totals are 301 leading, 354 trailing, and 330 symmetric. Their initial e, i, 
and (j) can be seen in Figs. |4| [s] and [6] Overall, librating Twotinos in the pre-run were 31.1% 
leading, 33.9% trailing, and 35.1% symmetric; Poisson statistics gives a precision of 0.8% for these 
proportions. There is a hint in Fig [6] that the phase space volume at low values of e < 0.1 
increases with inclination, indicating a change in the dynamical structure of the 2:1 resonance 
with inclination; it would be interesting to examine this in a future study with larger numbers of 
particles. 

A final piece of information obtained from the pre-runs is the amplitude of the short-period 
oscillations in a, e, and i caused by the three inner giant planets. These oscillations occur on 
timescales as small as tens of years, and are easily averaged over by the 2,000-yr output intervals 
of the pre-runs. Therefore, we define for each particle over each output interval a quantity da = 
{cLmax — 0-min)/'^ (and similarly de and di). We find that da, de and di remain remarkably constant 
over the 150-kyr pre-run integration length (see Fig. [s]). We calculated the mean values of da, de, 
and di over the 150-kyr interval for each particle, and adopt these as the amplitudes of the short 



period perturbations (see Section 3.2). 



The solid lines in Fig. |2]show the distribution of resonance libration amplitude, A(/>, from our 
pre-runs. These distributions obtain from a large number of particles sampling uniformly a wide 
range of e and i, and thus reflect the relative volume of resonance phase space characterized by 
any given value of A(j) — thus we may call it the volumetric distribution of Ac/). Note that this 
distribution includes all particles that are stable for 150,000 yr (a few libration periods), and thus 
does not reflect longer-term stability. From the shapes of these distribution functions, we can 
infer the relative volumes of phase space at various libration amplitudes. For Plutinos, we see 
that relative volumes within the resonance zone are nearly uniformly distributed in the interval 
45° < A(j) < 160°, as indicated by the nearly constant slope of the cumulative fraction, while 
values of A(j) < 45° represent less than ~ 10% of the resonance phase space. Leading and trailing 
Twotinos have most of their phase space volume in the range 20° < Ac/) < 75°; smaller and large 
libration amplitudes represent very little resonance phase space volume. For symmetrically-librating 
Twotinos, the phase space volume is mainly in the range 135° < Ac/) < 165°, though amplitudes as 
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Fig. 2. — The solid lines show the distribution of libration amplitudes in the preruns, an unbiased 
population which reflects relative volumes in resonance phase space. The dotted lines show the 
weighted distribution of libration amplitudes in the set of particles selected for our full 1-Gyr 
integrations. 



low as 90° are possible. 

With the initial conditions generated as described above, three numerical integrations were 
performed for a duration of 1 Gyr each. Run PO for Plutinos included only 3:2 resonant test 
particles in the gravitational field of the Sun and the 4 giant planets; similarly, the one Twotino 
run, TO, included only 2:1 resonant test particles. A second Plutino run, PI, included Pluto as an 
additional massive perturber for the 3:2 resonant test particles. 



3.2. Proper Elements 

The osculating orbital parameters of a resonant KBO can vary in a number of different ways, 
but the dominant variation is due to resonant perturbations of Neptune. The libration of the 



resonant argument (p is accompanied by correlated librations of a and e. In the approximation 
of the circular restricted 3-body model for resonant orbits (which pertains well to the qualitative 
characteristics of our particle orbits), the resonant librations of a and e are constrained by the 



conservation of the action A'^, also known as the "adiabatic invariant" (see, e.g., Murray and Dermott 



T999| ch. 8.12), 

(p + l)\/l-e2cosil , (2) 



N 



P 



where p = 1 for 2:1 resonant particles and p = 2 for 3:2 resonant particles. The librations of 
a and e are 90° out of phase with the librations of (p. N is conserved under resonant motions 
only; perturbations from massive bodies other than Neptune, and indeed non-resonant terms in 
the Disturbing Function due to Neptune, do not conserve A^. The resonant variations of these 
dynamical parameters have characteristic amplitudes and periods which depend on the particle's 
eccentricity and libration amplitude; the periods are in the range of 10^ to 10^ years. Any of these 
libration amplitudes (e.g., A4> or Aa) is diagnostic of how stable the particle's orbit is, with the 
smallest libration amplitudes usually being the most stable. However, numerical analysis of the 
resonant librations is not a simple matter, as other modes of variation occur simultaneously in each 
particle's evolution, including fast quasi-periodic variations as well as secular perturbations. 

Non-resonant perturbations from the 3 inner giant planets (Jupiter, Saturn, and Uranus) give 
rise to a suite of fast variations in a resonant particle's orbital elements. These variations can have 
significant amplitudes, but their periods are on the order of the period of conjunctions between 
the disturbing planet and the test particle. Non-resonant perturbations also arise from Neptune. 
For the Plutinos and Twotinos, these short period variations have periods ranging from ~ 12 yr 



(owing to Jupiter) to ~ TNep/[^ —p/{p + 1)] (Tiscareno 2004 1; the latter are owed to conjuctions 



with Neptune, and have periods 493 yr for the 3:2 resonance (p = 2) or 329 yr for the 2:1 resonance 
{p = 1). The amplitudes of these short-period variations remain very nearly constant and they 
average to zero over the resonance libration periods. 



The giant planets also give rise to long period secular variations in test particle orbits (Knezevic 
et al.|[l991 |. The semimajor axis is not affected by secular variations, but all other orbital elements 
are. In the outer solar system, these secular variations occur on ~Myr timescales, much slower than 
the mean-motion resonant librations. Thus the latter will continue to move along lines of constant 
"adiabatic invariant" N, even as the secular variations cause N itself to oscillate. 



Fig. ^ shows a sample particle from the pre-run (see Section 3.1 ) that demonstrates all of the 
behaviors described here. The resonant libration, with a period of 30,000 yr for this particle, is 
clearly visible in the time evolution of a, e, and (/>, while the evolution of i is dominated by the 
long-period secular variation. The secular variation is also perceptible as a gentle downward trend 
in the eccentricity evolution. The fast variations are indicated by the difference between the upper 
and lower solid lines (note that the output timestep for the pre-run is 2,000 yr, small enough that 
only the fast variations are smoothed over) , and the constancy of their amplitude is easily seen by 
the consonance of the two dotted lines in each panel. 
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Fig. 3. — Pre-run output for one test particle, illustrating the different kinds of variation of the 
orbital elements as discussed in Section [3 .21 Solid lines are recorded maximum and minimum values 
of a, e, and i over 2,000-yr intervals, and instantaneous values of the resonant argument. Dotted 
lines in the top three panels show the maxima and minima shifted by the constant values da, de, 
and di, and simply the midpoint value of A^ in the bottom panel. Note that this particle is a 
Twotino, librating in the trailing mode. 
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To tease out the diffusive evolution of particles due to long-term chaos, we adopt the method 



of Morbidelli (19971, as follows. In the full 1-Gyr integration, we record for each particle, at 1-Myr 
intervals, the minimum and maximum values of a, e, and i, over a sliding window with a length of 
10 Myr (thus the intervals are [0,10] Myr, [1,11] Myr, [2,12] Myr, etc.). Additionally, instantaneous 
values of 0,, uj, and M are also recorded every 1 Myr, so that we can follow directly the evolution 
of the resonance angles. Note that the results of this procedure will not explicitly display any 
of the behaviors described above and seen in Fig. |3] since all have periods smaller than 10 Myr; 
rather, these processes are smoothed over, and their amplitudes incorporated into the differences 
between the recorded maximum and minimum values. We then remove the amplitude from the fast 
oscillations induced by the giant planets, by subtracting from the maximum values (and adding 
to the minimum values) the short-period amplitudes da, de, and di that were obtained from the 
pre-runs (see Section |3.1[) . Finally, we remove secular oscillations from the proper eccentricities 



by readjusting them to correspond to the midpoint value of the adiabatic invarianlj^ A'^ over each 



10-Myr interval. A more detailed discussion of this so-called 'W algorithm" is given in Morbidelli 



(19971. 



The final proper elements defined by this process are the minimum (which we now refer to 
as ai, ei, and ii) and the maximum (a2, 62, and ^2) values that represent variation due only to 
the resonant librations. These would be stationary in the long term but for chaotic diffusion. The 
evolution of the the proper element set (ai, ei,ii) for four example particles is shown in Fig. [I] the 
chaotic evolution of these particles ranges from very stable (little change in the proper elements) 
to strongly chaotic (large and rapid changes in the proper elements). To quantify the chaotic 
diffusion, we determine the variation Aai between each pair of consecutive sliding-average time 
windows described in the previous paragraph. Then we calculate at each location in phase space 
the chaotic diffusion coefficient D = ((Aai)^)/Ai (see Section ^, where At = 10 Myr, with 
the average taken over all particles at all times at that phase-space location. We note that this 
calculation could equivalently have been done with 02. With this procedure, we find values of D 
in the range 10~^-10~^ AU^/Myr for the Plutino and Twotino particles in our 1 Gyr simulations. 

To calibrate the chaotic diffusion rates of the particles, we can compare them with the chaotic 
diffusion rates for the planets, which we expect to be relatively small. Numerical experiments have 



shown the evolution of the giant planets to be regular or only very weakly chaotic ( Hayes[2008 and 
references therein), and their orbital elements remain within narrow bounds over multi-gigayear 
timescales (|Laskar|1997[ |Ito and Tanikawa|2002 ) so that their proper elements ought to be very close 



to stationary. Table [T] shows the rms deviations of the proper elements for each of the giant planets 
over all of our 1-Gyr integrations. (Our values for the first three planets are in agreement with 



those quoted by Morbidelli (19971, while our values for Neptune are a factor of 3 to 5 larger. We 



do not know the source of the latter discrepancy, but the detailed study by Hayes (20081 suggests 



■""This "TV algorithm" is a 2-dimensional calculation performed in a — e space, which provides a good approximation of the 
secular variations in e. The complexity of adding a third dimension to this process precludes a similar calculation of the secular 
variations in i. However, we note that the latter are generally only a few degrees in amplitude. 
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Table 1: Rms deviations of proper orbital elements for the giant planets over all 1-Gyr integrations 
described in this paper 



Ua (AU) (Te Oi 

Jupiter 0.00003 0.00004 0.001° 

Saturn 0.00015 0.00020 0.004° 

Uranus 0.0014 0.00054 0.018° 

Neptune 0.0040 0.00032 0.006° 



that such differences are not surprising in practically similar long-term integrations of the outer 
solar system that may differ by adopting planets' initial conditions at different epochs and/or use 
different orbit integrators.) The magnitudes of the rms variations of the planets' elements provide a 
measure of the numerical resolution of our study for the detection of chaotic diffusion. The largest 
rms deviation, aa for Neptune, is equivalent to a diffusion coefficient D « 10"*^ AU^/Myr, which 
is three orders of magnitude smaller than the smallest diffusion coefficient we find for the resonant 
particles. 



4. Results 

Initial parameters and resonance residence times for particles in the three full integrations 
are shown in Figs. |4]~-[6j We define a Plutino as having escaped the resonance as soon as ai falls 
below 38.77 AU or 02 exceeds 40.17 AU; for Twotinos, the corresponding limits are 47.16 AU and 
48.56 AU. Both of these intervals are ±0.7 AU from the central resonant value; they are over- 
generous in that all resonant particles remain within the intervals, while some particles may remain 
briefly within the intervals despite being no longer in resonance. We find that once a particle has 
left the resonance for the first time, it spends very little time (on average, a few percent of each 
particle's dynamical lifetime) back in the resonant range of semimajor axis. Therefore, we consider 
a particle to be permanently non-resonant once it has first left the resonant range. 

Maps of the dynamical diffusion rate within the resonant region for each integration are shown 
in Fig.ul For these maps, we divided the proper element planes, (a, e) and (e, i), into small bins, and 
calculated the average value D of the diffusion coefficient at each phase-space location, averaging 
over all particles at all times whose initial proper elements lie in that bin at the beginning of a 
sliding-average time window (see Section ^ . The grey scale maps show how these average diffusion 
coefficient values vary across the resonance phase space. We have also defined the "characteristic 
diffusion time" associated with the chaotic diffusion, r = 1? jD^ where we use L = 0.2 AU as a 
characteristic lengthscale (half-width) for the resonant region; the scale bar in Fig. 17] shows the 
grey scaling for both D and r. 
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Finally, in Fig. |8j we show plots of the fraction of particles that survive in resonance as a 
function of time, and the distribution of the characteristic diffusion time. 

These results are discussed in detail below for the Plutinos and Twotinos separately. 
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Fig. 7. — Maps of mean dynamical diffusion coefficient D. Also noted on the scale bar is the 
characteristic time for diffusion associated with each diffusion coefficient, r = L'^/D, where we use 
L = 0.2 AU as a characteristic lengthscale (half-width) for the resonant region. The two small 
boxes in the PI plots show the extent of the variation in Pluto's minimum and maximum proper 
elements, respectively. 



4.1. Plutinos 

In Figs. |4] and |5] we see that long-term stable resonant orbits for Plutinos are prevalent at 
small lS.(j) s-iicl e values, and only slightly less common at higher inclinations. In the a-e planes 
(upper panels in Fig. u\, the lowest diffusion values are at a ~ ares and low e. In the e-i planes 
(lower panels in Fig. u|, low diffusion rate is found at low e for the entire range of inclinations, 
and the Kozai resonance — in which the argument of pericenter librates about 90° , thus enhancing 
stability (Kozai 19621 — is visible as a band of lower diffusion rate just to the right of the highest 



diffusion rates (0.2 < e < 0.3). 

By comparing the PI results with those for PO, the effects of Pluto can be assessed. Comparing 
Figs. |4] and [5] a number of particles become less stable when Pluto is present, but many also become 
more stable; Pluto's presence brings only a 3% net decrease in the mean particle lifetime of a particle 
in the resonance. Furthermore, very little difference can be discerned between the diffusion maps 
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Fig. 8. — (a) Fraction of particles remaining in resonance, as a function of time, for runs PO (solid), 
PI (dashed), and TO (dot-dash). Dotted lines show power-law fit to the last 0.5 Gyr of data, with 
extrapolation to 4 Gyr. (b) Relative frequency of characteristic diffusion times from maps in Fig.JTl 



for the two runs (Fig. Wh. However, one noticeable effect of Pluto is that the number of particles 
remaining in the Kozai resonance (identified by tracking the libration of uj) for the entire 1 Gyr is 
13% lower for run PI as compared with run PO. Because Pluto also occupies the Kozai resonance, 
Kozai librators will tend to have orbital elements similar to Pluto's. Thus, their encounters with 



Pluto produce greater perturbation because of the lower encounter velocities (Opik 19761. There 



is also a significant number of particles (10% of the total) that experience "Trojan" behavior with 
respect to Pluto at one time or another (i.e., either libration or slow circulation of A — Ap), though 
few particles (0.5% of the total) experience long-term Trojan libration. 

Pluto's overall effect on our test particle Plutinos is most clearly seen in the rate at which 
particles leave the resonance (Fig. ^). At the end of the 1-Gyr integration, 39% of the resonant 
particles in PO remain in resonance, compared with 37% in PI. The percentage projected to remain 
after 4 Gyr is also only slightly less for PI (27%) than for PO (28%), and the loss function exponents 
(Section 4.3) are similar (6pi = —0.554, while bpo = —0.556). 



4-2. Twotinos 

In Fig. [6] we see that long-term stability in run TO is correlated with smaller A<j) values; but 
here we find (in some contrast with the Plutinos) that moderate eccentricities (0.1 < e < 0.3) are 
more likely to be stable than low eccentricities, and that stability falls off visibly with increasing 
inclination. For inclinations greater than about 15°, we find no long term stable symmetric librators. 
These patterns are also seen in the diffusion maps (Fig. Uj, where we also note that the region of 
stability spreads over a somewhat wider interval in semimajor axis than is the case for Plutinos. 
An increase in stability associated with the Kozai resonance is discernible in the e-i plane (lower 
TO panel in Fig. u| as a faint band near e ~ 0.4, but it appears to be much less prominent in 
Twotinos than in Plutinos. 
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As shown in Fig. [Sp,, 24% of particles in run TO remain in resonance after 1 Gyr, about two- 



thirds the fraction found above for Plutinos, and the loss function exponent (Section 4.3) is steeper 
(^To = —0.768), yielding only a projected 15% remaining after 4 Gyr. This indicates that the 2:1 
resonance has less overall long-term stability than the 3:2. 



4-3. Resonant Population Decay 

Fig. |8^ shows, as a function of time, the fraction of particles in each run that have not yet 
escaped the resonance. For all runs, a downturn or "elbow" can be seen at about 100 Myr. This 



was also noted by Morbidelli (1997) in his study (which corresponds to our run PO), though the 



effect is more modest in our results. This "elbow" can be explained if the volume of phase space 
with characteristic diffusion times in the range 10® < r < 10^ yr were much larger than the volume 
with shorter diffusion times, so that particles would begin escaping from the former region only 
after 10® yr. This explanation is supported by the plot in Fig. [Sp, which shows the distribution of 
dynamical diffusion times: we see that ~ 80% of the resonance phase space volume has characteristic 
diffusion time in excess of 10® yr and shorted diffusion times are found in only 20% of the resonance 
phase space. 

We extrapolate the fraction, fres, of remaining resonant particles using a functional form 
oc t^, which is preferred by 



MorbidelH 



(1997) on theoretical grounds, although other functional 



/r, 

forms would also fit our data, such as fres oc logt (iHolman and Wisdom 19931). For each run, we 



have made a least-squares fit to the last 0.5 Gyr of our results to find the loss function exponent b, 
which we use to predict the fraction of particles remaining after 4 Gyr. For Plutinos, our best-fit 
exponent is bpo = bpi = —0.55; this is similar to [Morbidelli s estimate of 5 = —0.5. For Twotinos 
(which [Morbidelli did not investigate), our best-fit exponent is steeper, bxo = —0.77. 



The uncertainty in the loss function exponent b, and in the extrapolated fraction of particles 
remaining after 4 Gyr, can be estimated using standard error-propagation techniques. The standard 
deviation of the fitted data from the model (a linear fit in log-log space) is no more than 2% of the 
range covered by the data. The uncertainties for our quoted values of b are in the third significant 
figure. The uncertainties (calculated from the covariance matrix, and thus assuming that our choice 
of a linear model is correct) for our quoted fractions of particles remaining after 4 Gyr, extrapolated 
from our 1-Gyr simulations, are in the fourth significant figure, though we choose to quote only 
two significant figures. 



4.4- Resonance Escapees 

Particles that escape from resonance in our simulations enter other dynamical populations. 
The escaped particles spend roughly equal amounts of time as Centaurs and as scattered disk 
objects (SDOs). (The definitions of "Centaur" are somewhat variable in the literature; we adopted 
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the definition of Chiang et al. (2007) for these populations: objects that are neither Resonant nor 



Classical KBOs are Centaurs if their perihelion is interior to Neptune's orbit, and SDOs otherwise.) 
Escaped Plutinos spend somewhat more time (52% in run PI, 55% in run PO) as Centaurs, while 
escaped Twotinos spend somewhat more time (52%) as SDOs. Combined, escaped Plutinos and 
Twotinos spend less than 1% of their time as classical Kuiper Belt objects (CKBOs) — that is, with 
mean eccentricity (e) < 0.2 and mean Tisserand parameter with respect to Neptune (T/v) > 3. 
We find that ~ 40% of resonance escapees survive to the end of our 1-Gyr integration. These 
typically survive as SDOs, their stability often enhanced by "resonance sticking"; this may involve 



a long-term stay in a single resonance, or several resonances may be visited (Levison and Duncan 



19971. 



If a particle's perihelion migrates inside the orbit of Neptune while it is unprotected by a 



resonance, it then enters the population of Centaurs (Tiscareno and Malhotra 2003 Horner et al 



2004 di Sisto and Brunini 2007). Centaur dynamics are dominated by scattering due to the giant 



planets. Centaur orbital elements tend to diffuse nearly evenly throughout the planet-crossing 
region of parameter space. We find that 27% of particles that do not survive the 1-Gyr integration 
enter the inner solar system (that is, to heliocentric distance r < 5 AU); this is consistent with 



previous results (e.g., Tiscareno and Malhotra (2003)). We also see in our simulations that resonant 



KBOs make a significant contribution to the SDO population, as well as directly to the Centaurs. 



Discussion 



5.1. Comparison with Previous Work 



Morbidelli (1997) explored chaotic diffusion in the 3:2 resonance using methods very similar to 



ours. We observe all of the dynamical behaviors that he describes, including objects which slowly 
diffuse from the inner parts of the resonance phase space, escaping only after several billion years. 
We also reproduce fairly well his map of chaotic dynamical diffusion in the 3:2 resonance. We 
have expanded on his work by adding greater resolution, greater coverage of particles deep in the 
resonance, and also by exploring the influence of Pluto; we have also extended the study to the 2:1 
resonance. We confirm Morbidelli's observation of an increase in the rate of escape from resonance 
after t ~ 100 Myr (Fig. ^), and note that similar behavior appears in the 2:1 resonance as well. 
We ascribe this to a relatively large volume in phase space that has characteristic diffusion times 
on that order (Fig. Isl; 



Yu and Tremaine (1999) investigated a semi- analytical model based on a simplified resonance 



Hamiltonian, to explore the dynamics of Pluto and the Plutinos. All Plutinos in their model expe- 
rienced significant interactions with Pluto, but their initial conditions were restricted to particles 
with orbits very similar to Pluto's. In our more diverse sample, we do observe a set of behavioral 



classes that is similar to what Yu and Tremaine described, including persistent Trojan behavior 
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("tadpole" and "horseshoe" orbits, and orbits that transition between them), particles with a slow 
circulation of A — Xp, and others that experience only intermittent Pluto Trojan behavior. About 
24% of particles in our simulation experience at least intermittent commensurability with Pluto, 
but only 7% experience persistent circulation of A — Ap, and less than 1% are in tadpole or horseshoe 
orbits. 



Nesvorny and Roig (2000, 20011 investigated the stability of the 3:2 and 2:1 resonances, re- 
spectively. They used the digital filter method of obtaining proper elements, and they measured 
resonance stability using the Lyapunov exponent and the minimum distance from Neptune, rather 
than computing the dynamical diffusion coefficients as we have. Still, our maps of resonance sta- 
bility have broad agreement in the overall shape of the resonant region of phase space, with the 
most stable regions of the resonances centered on the resonant values of a, and an additional zone 
of stability associated with the Kozai resonance. 



Nesvorny et al. (2000) performed a short-term numerical integration of the known Plutinos 
in order to determine their present proper orbital elements, including the resonance libration am- 
plitude A(j). They found a surprising lack of observed Plutinos deep in the 3:2 resonance (i.e. at 
small A(j)), as well as at moderate eccentricities (i.e., near ep) in the Kozai resonance, which they 
attributed to the effects of Pluto. They also carried out numerical integrations of particles with 
initial eccentricities and/or inclinations very close to those of Pluto, showing that Pluto ejects a 
large fraction which would otherwise be stable. We find a much smaller influence of Pluto upon 



the Plutinos, likely due primarily to [Nesvorny et al. s reliance on low- amplitude librators — all of 



their particles begin with Acp < 20°, and most with A(p < 10° — which, as they demonstrate, are 
much more likely to interact with Pluto. By contrast, our pre-runs (which do not include Pluto) 
indicate that less than 2% of phase space is represented by such low- amplitude librators, a fraction 
reflected in the initial conditions of our models. Consequently, we find Pluto's overall effects to be 
on the order of a few percent increase in objects lost from the resonance over 4 Gyr, rather than 
50%. This large difference is mainly because we adopt an initial libration amplitude distribution 
which is proportional to the phase space volume at different libration amplitudes. 



Our present work is a condensed but updated version of the study by Tiscareno ( 2004 ) , who 
additionally investigated the effects of a massive perturber embedded in the 2:1 mean- motion 
resonance. The earlier work contained an error in the planetary initial conditions and also used a 



less robust process for selecting initial conditions for the test particles. Tiscareno (2004 1 found a 
reduced impact of Pluto upon the Plutinos compared to previous work, but here we find Pluto's role 
even further reduced, likely due primarily to truly random selection of particles from the pre-run 
to avoid bias towards low-Ai?i librators, and secondarily to randomly distributed initial values of to 
and O rather than constraining them to equal Pluto's initial values. 
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5.2. Population Estimates 



Kavelaars et al. (20081 estimated that Plutinos are currently ~ 20% as abundant as classical 



Kuiper Belt objects. Kavelaars et al. do not give an estimate for the current population of Twotinos, 



but Chiang and Jordan (20021 estimated that Plutinos are ~ 3 times as abundant as Twotinos. 



Using our results from Sections 4.1 and 4.2 that 37% of the 3:2 resonance population survives 
for 1 Gyr (run PI) and 24% survives in the 2:1 resonance (run TO), and extrapolating the loss 
function curves to 4 Gyr (Section |4.3[ ) , we estimate that the Plutino population has decayed ~ 73% 
over the past 4 Gyr and the Twotino population has decayed ~ 85% over that time period. (The 
uncertainties in these estimates are small - see Fig. [8^ and discussion in Section |4.3[ future work 
could check these estimates by means of full 4-Gyr simulations.) In contrast with their present 
abundance, the Plutinos may have formerly been as much as ~ 75% as abundant as CKBOs, and 
Twotinos as much as ~ 45% as abundant. Thus the resonant population comprised a much larger 
fraction of the Kuiper Belt 4 Gyr ago, likely more than half of the total (classical + resonant). 
Furthermore, the Plutino/Twotino ratio 4 Gyr ago was about half what it is today. These relative 
population estimates for 4 Gyr ago are consistent with models that predict the two populations 



to have comparable numbers at the conclusion of Neptune's migration (e.g., Malhotra 1995 Hahn 
and Malhotr"a||2005[ )f] 



The above estimates and conclusions are based on the assumption of an initially volumetric 
distribution of A0, meaning that the frequency of a particular value of A(^ is proportional to the 
volume in phase space characterized by that value (particles that do not survive the first 1 Myr 
of the full integration are not counted, as they are not likely to have survived any process by 
which the Kuiper Belt could have been formed). It is unclear whether this assumption is robust; 
current alternative models of Kuiper belt structure formation have not noted the resonance libration 
amplitude distributions. If, instead of the volumetric distribution, only the very stable particles 
survive the formation process to become the initial resonant populations, then particle escapes 
overall would be much less common; consequently, the fraction of resonant particles that survive 
to the present age of the solar system would be much larger, and neither of the two conclusions 
above would be very robust. For illustration, we consider the relatively extreme case that initially 
only the small amplitude librating orbits, Ac/) < 30°, were populated; in this case, we estimate 
from our simulations that approximately 67%, 56% and 38% of the populations in our model 
runs PO, PI and TO, respectively, would survive to 4 Gyr (instead of approximately 28%, 27%, 
and 15%, as described above). In this case, the projected primordial resonant populations and 
the Plutino/Twotino ratio are significantly different, and, as discussed above, Pluto's perturbing 
effects are significantly larger. These considerations show that models of Kuiper Belt structure 
formation should pay attention to the resonance libration amplitude distributions, as the latter are 
critical in determining the history of the resonant populations; conversely, the resonance libration 



^It is possible that the early Kuiper Belt history may have included rapid loss from the strongly unstable regions of the 
resonances, and such loss is unconstrained by the presently observed populations. 
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amplitude distributions may be critically useful in constraining alternative theoretical models. This 
is a possible direction for future research. 



6. Conclusions 

We carried out extensive numerical integrations to probe the characteristics of the 3:2 (Pluti- 
nos) and 2:1 (Twotinos) mean motion resonances with Neptune. In each case we obtained a 
representative uniform sample of initial conditions in the resonance phase space, and integrated 
these for 1 Gyr. Our results can be summarized as follows: 

1. In both resonances, for orbits that are stable on resonance libration timescales (~ 10^ yr), 
the resonance phase space volume is described by the following distributions of the libration 
amplitudes (A^). For Plutinos, /S.(j) has uniform distribution between 45° and 160°; only a 
small fraction, < 5%, of the phase space volume has At/) < 30°. For Twotinos, symmetric 
librators have a nearly uniform distribution of Ac/) between 145° and 165°, while asymmetric 
librators (both leading and trailing) have a nearly uniform Ac^ distribution between 20° and 
75°; smaller and larger libration amplitudes represent < 10%. 

2. The presence of Pluto has only modest effects on the Plutino population. It disrupts the 
stability of objects in the Kozai resonance and causes some particles to undergo libration in 
a 1:1 ("Trojan") resonance. The overall fraction of Plutinos that survive for 1 Gyr is 39% 
when Pluto's perturbations are neglected, and 37% when they are included. Projected to 
4 Gyr, these two surviving fractions are 28% and 27%, respectively (Fig. [8^). If, however, 
the initial distribution of A</) were heavily weighted towards low-amplitude librators, Pluto's 
effects would likely be more important. 

3. The 2:1 resonance has weaker overall long-term stability than the 3:2; only ~ 24% of Twotinos 
survive in our 1 Gyr integration, and only ~ 15% are projected to survive for 4 Gyr (Fig.^ 



4. The population decay rates obtained in our models indicate that, 4 Gyr ago, resonant KBOs 
made up more than half of the Kuiper Belt, and the Plutino/Twotino ratio was closer to unity 
than it is now. These estimates depend on our assumption of an initial distribution of particles 
proportional to the local phase space volume in the resonance and to the extrapolation from 
our 1-Gyr integrations. These population ratio estimates for the 4-Gyr old Kuiper belt are 
roughly consistent with simple models of Neptune's migration and resonance sweeping of 



the Kuiper Belt (Malhotra 1995), and provide a useful constraint for more detailed models 



(Morbidelh and Brown 2004 Morbidelh et al. 20081. 
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